home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 20 / Cream of the Crop 20 (Terry Blount) (1996).iso / math / alged34.zip / ALGEDSRC.ZIP / ALGCALC.C next >
C/C++ Source or Header  |  1996-06-06  |  17KB  |  556 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. /*-----------------------------------------------------------------
  11.    prime factor - replace a number with prime factors
  12. */
  13. void primefact(node *p) {
  14.   int i;
  15.   double v,f;
  16.  
  17.   for (i=0; i<p->nump; ++i)
  18.     primefact(p->parm[i]);
  19.  
  20.   if (p->kind==NUM && !!(v=fabs(p->value)) && whole(v)) {
  21.     if (p->value<0) {
  22.       p->kind = MUL;
  23.       p->nump = 2;
  24.       p->lf = newnum(v);
  25.       p->rt = newnum(-1);
  26.       p = p->lf;
  27.     }
  28.     for (f=2; f<=maxrat && f<v; ++f)
  29.       if (whole(v/f)) {
  30.         v /= f;
  31.         p->kind = MUL;
  32.         p->nump = 2;
  33.         p->lf = newnum(v);
  34.         p->rt = newnum(f);
  35.         p = p->lf;
  36.         f = 1;         /* start loop at the beginning */
  37.       }
  38.   }
  39. }
  40.  
  41. /*--------------------------------------------------------------------
  42.    rational search - this searches for a denominator to a number.
  43. */
  44. int rational_search(node *p) {
  45.   int i,r=0,s=1;
  46.   static double v,n,x,d,n1,d1,x1;
  47.   double err = pow(10,-sigdig);
  48.  
  49.   twirl();
  50.   if (p->kind==NUM) {
  51.     v = fabs(p->value);
  52.     if (p->value<0) s=-1;        /* sign */
  53.     x = modf(v,&n);
  54.     if (x > err) {              /* v is not already an integer */
  55.       x1=5;
  56.       for (d=2; d<=maxrat; ++d) {
  57.         modf(d*v+0.5,&n);            /* n is the round(d*v) */
  58.         x = fabs(n/d - v);           /* x is the error */
  59.         if (x < x1) {
  60.           n1 = n; d1 = d; x1 = x;
  61.         }
  62.       }
  63.       if (x1 < err) {          /* convert p to a ratio of n/d */
  64.         if (n1==0) {
  65.           p->value = 0;
  66.         }
  67.         else {
  68.           p->kind = DIV;
  69.           p->nump = 2;
  70.           p->lf = newnum(n1*s);
  71.           p->rt = newnum(d1);
  72.         }
  73.         ++r;
  74.       }
  75.     }
  76.     else {                /* throw away the very small part */
  77.       p->value = n*s;
  78.       ++r;
  79.     }
  80.   }
  81.   return r;
  82. }
  83.  
  84.  
  85.  
  86. /*-----------------------------------------------------------------
  87.    reduce a/b to lowest terms (with positive denominator)
  88. */
  89. int reduce(double *a,double *b) {
  90.   int r=0;
  91.   double i;
  92.  
  93.   if (*b<0) { *a=-*a; *b=-*b; }
  94.   for (i=2; i<=maxrat && i<=fabs(*a) && i<=*b; ++i)
  95.     if (whole(*a/i) && whole(*b/i)) {
  96.       *a /= i; *b /= i;
  97.       ++r; i=1;
  98.     }
  99.   return r;
  100. }
  101.  
  102. /*---------------------------------------------------------------------------
  103.    ration - this converts all the numbers to ratios of integers if
  104.             possible.
  105. */
  106. int ration(node *p) {
  107.   static double v,u,h,n1,d1,w9,pf,pp,rp;
  108.   static int m,j,k,w,n,f;
  109.   int i,r=0;
  110.   char s[30];
  111.  
  112.   for (i=0; i<p->nump; ++i)
  113.     r+=ration(p->parm[i]);
  114.  
  115.   if (p->kind==NUM && !whole(p->value)) {
  116.     u = fabs(p->value);
  117.     v = frexp(u,&m);             /* u ==> v * 2**m */
  118.     if (m > 0) {
  119.       n = sigdig - m*log10(2);   /* convert m to base 10 logarithm */
  120.       if (n < 1) {
  121.         modf(p->value,&p->value);  /* truncate to integer */
  122.         return r+1;
  123.       }
  124.       v = modf(u,&h);
  125.       m = 0;
  126.     }
  127.     else {                   /* u is a small number */
  128.       h = 0;
  129.       v = modf(u,&h); m=0;  /* suppress tiny fractions */
  130.       n = sigdig;
  131.     }
  132.     sprintf(s,"%1.18f",v);      /* we know 0 < v < 1 */
  133.     /* Look for repeating pattern in s */
  134.     for (f=0; f<n; ++f) {
  135.       k = n-f-1;
  136.       for (w=1; w<k; ++w) {
  137.         for (i=0; i<w; ++i) {
  138.           for (j=f+i+w; j<n && s[j+2]==s[f+i+2]; ) j+=w;
  139.           if (j<n) break;    /* failed */
  140.         }
  141.         if (i==w) break;   /* success */
  142.       }
  143.       if (w<k) break;   /* success */
  144.     }
  145.     if (w<k) {          /* success */
  146.       s[f+w+2] = 0;
  147.       sscanf(s+f+2,"%lf",&rp);     /* repeating part */
  148.       s[f+2] = 0;
  149.       s[1] = '0';
  150.       sscanf(s+1,"%lf",&pp);       /* prefix part */
  151.       w9 = pow(10,w)-1;           /* w nines */
  152.       pf = pow(10,f);
  153.       n1 = (pf*h + pp)*w9 + rp;
  154.       d1 = w9*pf*pow(2,-m);
  155.       if (p->value<0) n1 = -n1;
  156.       reduce(&n1,&d1);
  157.       if (n1==0) p->value = 0;
  158.       else {
  159.         p->kind = DIV;
  160.         p->nump = 2;
  161.         p->lf = newnum(n1);
  162.         p->rt = newnum(d1);
  163.       }
  164.       ++r;
  165.     }
  166.     else             /* try to convert using the search method */
  167.       r += rational_search(p);
  168.   }
  169.   return r;
  170. }
  171.  
  172.  
  173. /*--------------------------------------------------------------------
  174.    Move numbers within clag.
  175. */
  176. int movenums(node *p,int up,int oper) {
  177.   int i,r=0;
  178.   node *a,*b;
  179.  
  180.   for (i=0; i<p->nump; ++i)
  181.     r+=movenums(p->parm[i],up,oper);
  182.   if (p->kind!=oper) return r;
  183.   a = p->lf;
  184.   b = p->rt;
  185.   i = 1;
  186.   if (a->kind!=oper) {
  187.     a=p; i=0;
  188.   }
  189.  
  190.   if (a->parm[i]->kind==NUM) {
  191.     if (p->rt->kind==NUM) {           /* combine nums */
  192.       if (oper==ADD) a->parm[i]->value += p->rt->value;
  193.       else           a->parm[i]->value *= p->rt->value;
  194.       a = p->lf;
  195.       nodecpy(p,a);
  196.       freenode(a);
  197.       freenode(b);
  198.       ++r;
  199.     }
  200.     else if (up) {                         /* switch */
  201.       p->rt = a->parm[i];
  202.       a->parm[i] = b;
  203.       ++r;
  204.     }
  205.   }
  206.   else if (b->kind==NUM && !up) {   /* switch */
  207.     p->rt = a->parm[i];
  208.     a->parm[i] = b;
  209.     ++r;
  210.   }
  211.   return r;
  212. }
  213. /*-----------------------------------------------------------------
  214.    is complex?  returns >0 and a and b if p is a complex (1=real)
  215.    a is the imaginary part.  sorry, this is against convention!
  216. */
  217. int iscomplex(node *p,double *a,double *b) {
  218.  
  219.   if (aop(p->kind) && p->rt->kind==NUM &&              /* ai+b */
  220.       p->lf->kind==MUL && p->lf->rt->kind==VAR &&
  221.       p->lf->lf->kind==NUM && !strcmp(p->lf->rt->name,"i")) {
  222.     *a = p->lf->lf->value;
  223.     *b = p->rt->value;
  224.     if (p->kind==SUB) *b = -*b;
  225.     return 5;
  226.   }
  227.   if (p->kind==MUL && p->rt->kind==VAR &&              /* ai */
  228.       p->lf->kind==NUM && !strcmp(p->rt->name,"i")) {
  229.     *a = p->lf->value;
  230.     *b = 0;
  231.     return 3;
  232.   }
  233.   if (aop(p->kind) && p->rt->kind==NUM &&              /* ia+b */
  234.       p->lf->kind==MUL && p->lf->lf->kind==VAR &&
  235.       p->lf->rt->kind==NUM && !strcmp(p->lf->lf->name,"i")) {
  236.     *a = p->lf->rt->value;
  237.     *b = p->rt->value;
  238.     if (p->kind==SUB) *b = -*b;
  239.     return 5;
  240.   }
  241.   if (p->kind==MUL && p->lf->kind==VAR &&              /* ia */
  242.       p->rt->kind==NUM && !strcmp(p->lf->name,"i")) {
  243.     *a = p->rt->value;
  244.     *b = 0;
  245.     return 3;
  246.   }
  247.   if (aop(p->kind) && p->lf->kind==VAR &&            /* i+b */
  248.       p->rt->kind==NUM && !strcmp(p->lf->name,"i")) {
  249.     *a = 1;
  250.     *b = p->rt->value;
  251.     if (p->kind==SUB) *b = -*b;
  252.     return 4;
  253.   }
  254.   if (p->kind==VAR && !strcmp(p->name,"i")) {          /* i */
  255.     *a = 1;
  256.     *b = 0;
  257.     return 2;
  258.   }
  259.   if (p->kind==NUM) {                              /* b */
  260.     *a = 0;
  261.     *b = p->value;
  262.     return 1;
  263.   }
  264.   return 0;
  265. }
  266.  
  267. /*-----------------------------------------------------------------
  268.    complex base raised to a power - convert to r^x * cos tx + i sin tx
  269.    b is the imaginary part.  p->rt is REUSED and p->lf is freed.
  270.     -----------------------------------------------------------------
  271.          Note: I am not using this function because (1) it is
  272.          questionable whether the cos+isin expansion is "simpler" than
  273.          exp(i*x), and (2) because it was happening needlessly when the
  274.          complex base was not completely sorted and simplified.
  275.  
  276. void complexpow(node*p,double b,double a) {
  277.   node *c,*s;
  278.   double t;
  279.  
  280.   t = atan2(b,a);         //  -pi < t <= pi
  281.   c = newnode();
  282.   c->kind=FUN;
  283.   strcpy(c->name,"cos");
  284.   c->nump=1;
  285.   c->lf = cons(deepcopy(p->rt),MUL,newnum(t));
  286.  
  287.   s = newnode();
  288.   s->kind=FUN;
  289.   strcpy(s->name,"sin");
  290.   s->nump=1;
  291.   s->lf = deepcopy(c->lf);
  292.  
  293.   freetree(p->lf);
  294.   p->kind = MUL;
  295.   p->lf = cons(newnum(hypot(a,b)),EXP,p->rt);
  296.   p->rt = cons(c,ADD,con